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NOMENCLATURE 

p   =  density,  Kg/m3 

u  =  one  dimensional  velocity,  m/sec 

w  =  crack  opening  displacement,  m 

v  =  gas  seepage  velocity  normal  to  crack,  m/sec 

e  =  specific  energy,  J/m3 

C  =  constant  volume  specific  heat,  J/Kg-K 

m  -  mass  flow  rate,  Kg/sec 

P  =  gas  pressure,  Pa 

q"=  heat  transfer  into  crack  wall,  W/m2 

T  -  temperature ,  K 

f  =  fanning  friction  factor 

a   =  stress, Pa 

k  =  permeability,  m2 

$  =  porosity 

/x  =  viscosity,  kg/m-sec 

t  =  time,  sec 

Pr=  Prandtl  number 

q  =  thermal  diffusivity,  m2/sec 

K  =  thermal  conductivity,  J/m-K-sec 

€  =  fracture  roughness,  m 
Re=  Reynolds  number 
n  =  burn  rate  exponent 


NOMENCLATURE 


B  =  burn  rate  constant,  m/sec-Pa 
D  =  propellant  web  thickness,  m 
C  =  propellant  covolume ,  m3/Kg 
F  =  propellant  energy  density,  J/Kg 

KT  TT=  stress  intensity  factors,  ?ajm 
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CHAPTER  1 
INTRODUCTION 

Oil  and  gas  wells  have  been  stimulated  using  explosives  since 
the  late  1800' s.  Early  efforts  tended  to  emphasize  the  potential  of 
high  explosives  placed  in  the  wellbore  to  create  an  almost 
instantaneous  source  of  high  pressure  gas  resulting  in  fractures 
around  the  borehole. 

Unfortunately,  using  explosives  may  result  in  creating  an 
intensely  damaged  zone  near  the  borehole  that  is  only  rarely 
associated  with  the  long,  radial  fractures  required  for  successful 
stimulation.  A  simple  explanation  is  when  detonation  occurs,  nearby 
rock  may  yield  and  plastically  deform.  When  the  stress  wave  passes, 
the  rock  unloads  elastically.  This  leaves  an  increased  borehole 
diameter  and  a  residual  stress  field,  sometimes  called  a  stress  cage, 
which  is  compressive  near  the  borehole.  The  existence  of  residual 
stress  regions  around  boreholes  that  have  been  subjected  to  explosive 
detonations  is  a  well -documented  phenomenon  (Schmidt,  et  al . ,  1981). 
Some  of  these  observations  have  been  made  as  a  result  of  field 
experiments  at  the  Nevada  Test  Site  by  Sandia  National  Laboratories. 

Multiple  fractures  are  highly  desirable  in  formations  which  have 
naturally  occuring  fractures  such  as  Devonian  shale.  This  can  be  seen 
in  Figure  1.1.  The  production  from  the  unstimulated  well  depends 
primarily  on  the  number  of  fractures  intersected  during  drilling. 
Hydraulic  fracturing  typically  produces  a  single  fracture  that  is 
likely  to  run  parallel  to  most  of  the  existing  fractures  since  its 
orientation   is  governed  by  in-situ  stresses  that  probably  also  govern 
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the  pattern  of  the  natural  fractures.  The  advantage  of  multiple 
fractures  is  that  although  they  may  not  extend  as  far  as  a  hydraulic 
fracture,  they  may  connect  the  well  to  naturally  occuring  fractures. 

An  alternative  to  explosives  is  to  tailor  the  pressure- time 
behavior  of  the  explosive  or  a  suitable  propellant  so  as  to  keep  the 
peak  pressure  and  the  loading  rate  below  the  level  that  would  result 
in  a  stress  cage.  Unfortunately,  the  combination  of  loading 
parameters  that  will  produce  optimal  multiple  fracturing  and  avoid  the 
formation  of  a  stress  cage  is  not  well  known. 

When  compared  to  conventional  hydraulic  fracturing,  tailored- 
pulse  loading  has  several  possible  advantages.  Hydraulic  fractures, 
which  are  propagated  at  pressures  that  are  slightly  higher  than  the 
minimum  in- situ  stress  and  pumping  times  that  are  on  the  order  of 
hundreds  of  seconds  to  hours,  typically  produce  only  a  single  fracture 
whose  orientation  is  perpendicular  to  the  miminum  principal  stress. 
The  higher  wellbore  pressures,  which  are  achieved  in  tailored-pulse 
loading,  are  needed  to  drive  cracks  in  less  favorable  directions  with 
respect  to  the  in-situ  stresses  in  order  to  produce  multiple 
fractures . 

In  addition,  tailored-pulse  loading  with  propellants  is  likely 
to  produce  little  formation  damage  due  to  the  interaction  of  the 
working  fluids  with  the  rock.  Little  water  is  produced  by  these 
propellant  materials,  and  the  products  have  very  little  time  to  react 
with  the  rock.  Some  hydraulic  fracturing  fluids,  however,  can  cause 
swelling  of  the  rock.  Economic  considerations  may  also  make  the 
techniques  more  feasible  since  igniting  an  explosive  or  a  propellant 
charge  in  a  well  can  be  relatively  inexpensive.   Very  little  equipment 


or   time   may  be   required  when  compared   to  even  a  small  hydraulic 
fracture  job. 

If  tailored  pulse  loading  is  to  be  successful,  the  fractures 
generated  must  extend  a  significant  distance  from  the  wellbore.  It 
has  been  estimated  (  Jean-Claude  Rogiers ,  1987  )  that  lengths  on  the 
order  of  100-200  meters  are  required.  Alternatively,  if  the  procedure 
is  used  only  as  a  clean  up  tool  for  damaged  boreholes,  much  shorter 
lengths  would  be  acceptable. 

1 . 1  Experimental  Data 

Experiments  relating  to  the  multi- fracturing  process  have  either 
been  full  field  experiments  or  scaled  down  laboratory  tests.  At  the 
Sandia  National  Laboratories,  Schmidt  et  al.  (1981)  performed  a  series 
of  five  full-scale  tests  to  evaluate  various  multi-fracture  concepts. 
The  tests  were  conducted  at  the  Nevada  Test  Site  in  horizontal 
boreholes  drilled  from  a  tunnel  in  ash- fall  tuff. 

Various  propellants  and  explosives  were  tested.  After  the 
shots,  the  rock  was  mined  back  along  the  borehole  which  allowed  direct 
observation  of  the  fracture  patterns.  Sandia' s  results  led  to  a 
better  understanding  of  the  pressure  loading  rates  required  for 
successful  stimulation.  It  should  be  noted  that  no  field  tests  have 
observed  fractures  greater  than  about  10  meters. 

Previous  laboratory  work  on  multiple  radial  fracturing  is  not 
extensive,  because  of  the  difficulties  of  working  with  realistic 
energy  sources  on  a  lab  scale  and  making  appropriate  pressure 
measurements  and  fracture  observations.  Through  laboratory  scale 
tests   it  has  been  shown  that  multiple  fractures  do  occur  under  proper 


conditions.  This  can  be  seen  in  Swift  et  al .  (1981),  and  Fourney  et 
al.  (1981,  1983).  However,  small  scale  experiments  are  not  adequate 
to  address  the  fundamental  problem  of  crack  length. 

1 . 2  Existing  Analysis 

The  most  complete  gas  model  that  has  been  developed  with  the 
gas -driven  fracture  problem  in  mind  is  the  work  of  Nilson  and 
Griffiths  [Nilson,  et  al . ,  (1985,  1986),  Griffiths,  et  al . ,  (1986)]. 
In  their  approach,  similarity  solutions  are  obtained  for  the  problem 
of  a  planar  gas-driven  fracture  propagating  into  a  brittle  elastic 
solid.  Local  values  of  the  fluid  pressure,  temperature  and  speed  are 
computed  from  the  one-dimensional  mass,  momentum  and  energy  equation 
governing  the  turbulent  flow  of  a  ideal  gas.  The  equations  used  are 
very  similar  to  the  ones  that  will  be  derived  in  Chapter  2.  Crack 
opening  displacements  and  the  stress  intensity  at  the  crack  tip  are 
calculated  from  the  gas  pressure  distribution,  in  accordance  the 
quasi- steady  integral  relations  from  linear/elastic  fracture 
mechanics . 

The  geometry  is  presumed  to  be  either  planar  or  axisymmetric . 
The  appropriate  assumption  depends  on  the  relative  length  of  the 
fractures  compared  to  the  pressurized  length  of  the  borehole.  Opening 
displacements  are  calculated  from  the  theory  of  linear  elasticity 
(Sneddon  and  Lowengrub ,  1969)  using  the  following  quasi -steady 
formula: 

ww  nG        Jx  J0  2      2.1/2    ^R+£;CH  (.2      2    1/2 


where  P  is  the  internal  pressure  which  varies  along  the  fracture,  R  is 
the  radius  of  the  cavity,  G  and  u  are  the  shear  modulus  and  Poisson's 
ratio,  and  a  is  the  in-situ  stress  acting  normal  to  the  plane  of  the 
fracture.  The  aperture,  w,  at  any  given  location,  x,  depends  upon  the 
entire  pressure  distribution  along  the  fracture. 

The  stress  intensity  at  the  tip  of  the  fracture  also  depends 
upon  the  pressure  loading  everywhere  along  the  fracture.  Propagation 
occurs  when  the  stress  intensity  factor  reaches  the  critical  value. 
The  propagation  may  be  viewed  as  a  equilibrium  process,  in  which  the 
stress  intensity  is  maintained  at  the  critical  value  as  the  fracture 
propagates . 

Nilson,  in  the  above  described  work,  leaves  out  inertia  terms  in 
the  conservation  of  momentum  equation.  Swenson  and  Taylor  (1983) 
developed  a  finite  element  code  that  uses  a  smeared  crack  approach. 
This  model  doesn't  have  the  capability  to  model  gas  flow  in  the 
fracture.  CRACKER  (Swenson,  1985)  allows  the  crack  to  be  a  distinct 
part  of  the  mesh,  with  finite  elements  on  both  sides  of  the  crack 
face.  This  allows  a  clean  implementation  of  fracture  mechanics 
concepts.  However,  the  pressure  distribution  along  the  crack  face  is 
not  known.   Therefore  a  pressure  profile  is  assumed. 

1 . 3  Objectives  and  Scope 

The  primary  objective  of  this  thesis  is  to  develop  a  fully 
coupled  dynamic  model  of  gas -driven  fractures  and  then  to  use  this 
model  to  address  the  fundamental  question  of  maximum  obtainable  crack 
lengths.  This  will  allow  a  more  realistic  evaluation  of  tailored 
pulse   loading.   In  addition,  the  model  of  gas-driven  dynamic  fracture 


provides  a  general  tool  to  be  used  in  other  areas  such  as  rock 
blasting  and  containment  of  nuclear  blasts.  This  objective  was  met  by 
coupling  a  model  of  compressible  gas  flow  to  an  existing  finite 
element  code  which  models  dynamic  fracture  (Swenson,  1985)  .  The 
program,  including  the  modifications  made  for  this  thesis,  is 
completely  menu  driven  with  preprocessing,  analysis,  and 
postprocessing  available  to  the  user. 

The  remainder  of  this  thesis  is  arranged  in  five  chapters. 
Chapter  2  shows  the  derivation  of  all  equations  and  the  implementation 
into  finite  difference  form.  Chapter  3  is  a  discussion  on  the 
subroutine  GAS_SOLVE  and  some  of  the  interactive  features  available. 
Chapter  4  briefly  describes  the  problems  used  to  verify  the  model. 
Chapter  5  outlines  the  results  found  and  Chapter  6  gives  some 
conclusions  and  recommendations  for  future  work. 


CHAPTER  2 


FINITE  DIFFERENCE  GAS  DYNAMICS 


2 . 1  Governing  Equations 

Consider  the  elemental  control  volume  shown  in  Figure  2.1.  We 
can  derive  the  continuity  equation  by  doing  a  mass  balance  on  this 
element.  For  this  analysis  we  will  assume  compressible,  one 
dimensional  unsteady  flow.  There  will  be  no  property  variations  in 
the  direction  perpendicular  to  flow.  The  mass  entering  the  left  face 
of  the  element  per  unit  time  is: 

puw  (2.1) 

where  p  is  the  density,  u  is  the  one  dimensional  velocity,  and  w  is 
the  crack  opening  displacement.  The  mass  leaving  the  right  face  of 
the  element  per  unit  time  is: 

8 (gUW) ,  /o  o\ 

puw  +  — r* — 'dx  (2.2) 

dx  v    ' 

The  mass  leaving  the  top  and  bottom  faces  of  the  element  per  unit  time 
is : 

2pvdx  (2.3) 

where  v  is  the  seepage  velocity. 
The  time  rate  of  change  of  mass  inside  the  element  is: 

^T  dx  -  (w  at  +  pdi  )dx  (2-4) 

where  ~r   is  the  crack  opening  velocity.   A  mass  balance  on  this  system 


yields  the  conservation  of  mass  equation: 

ifi  .    1 

dt    "  w 


■  (  *£»«)  +  2pv  +  pft    }         (2.5) 


puw 


dx 


Figure  2.1:  Differential  Element  for  Conservation  of  Mass 
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Figure  2.2:  Differential  Element  for  Conservation  of  Energy 


Using   the   elemental   control   volume   shown   in  Figure  2.2  the 
conservation   of   energy   equation   is   derived.    We  will   assume 
constant  viscosity,   thermal  conductivity,  and  specific  heats.   Then, 
for   the   element   shown,    the   energy  balance  may  be  written: 
Energy  convected  in  left  face  -  energy  convected  out  the 
top  and  bottom  faces  +  net  work  done  on  the  element 
-  energy  convected  out  the  right  face  =  time  rate  of  change 
of  energy  in  the  element 
The  energy  entering  the  left  face  per  unit  time  is: 

epuw  (2.6) 

u2 

where   e=CT+x   ,  T  is  the  temperature,  and  C   is  the  constant 
v      2   '  ^        '       v 

volume   specific  heat.    The   flow  work  done  on  the  left  face  can  be 
found  in  the  following  manner: 

The  volume  flow  rate  (Q)  =  — 

P 

The  velocity  (u)  =  7  =  — 

J  A    pw 

The  force  applied  on  the  left  face  =  Pw 

The   flow  work  per  unit   time   (W)  is  the  product  of  the  force 
exerted  and  the  velocity: 

W  =  2!  (Pw)  =  m  p  (2  7) 

pwv      p  N    ' 

The   energy   leaving   the   right   face   of  the  element  per  unit  time: 

d (epuw)  ,  ,_  .. 

epuw  +   ' »K — xdx  (2.8) 
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The  flow  work  done  on  the  right  face  per  unit  time  is: 

A£  +  f(lhf)  (2.9) 

Energy  leaving  the  top  and  bottom  faces,  due  to  seepage,  per  unit  time 
is : 

2epvdx  (2.10) 

Time  rate  of  change  of  energy  inside  the  element: 

8 ( Ew )  ,     ,   5E   r3w  XJ 

IT  dx  =  (w  at  +  Eat  )dx 

where  E  =  ep .   The  work  done  to  expand  the  crack  and  by  friction: 

W  =  -P  fr  dx  +  2rAu  (2.11) 

dt 

Writing   the   energy  balance  corresponding  to  the  quantities  shown  in 

Figure  2.3  yields  the  conservation  of  energy  equation: 

1    3(uw(E  +  P)  +    +  p   +  2v   +  p)+2(^„  .  2rAu  j  _  dg      (2.12) 
w        dx  \     /   i  at 

where  q"  is  the  heat  transfer  into  the  crack  wall. 

Doing   the   same   type  of  analysis  on  the  element  shown  in  Figure 
2 . 3  we  can  find  an  expression  for  the  conservation  of  momentum. 
The  force  balance  on  the  element  is: 

SF  =  increase  in  momentum  flux  +  rate  of  change  of 
momentum  in  the  element 
The   momentum   flux  in  the  x  direction  is  the  product  of  the  mass  flow 
through  a  particular  side  of  the  control  volume  and  the  x  component  of 
velocity   at   that  point.    The   mass   entering   the  left  face  of  the 
element  per  unit  time  is: 

puw  (2.13) 
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Figure  2.3:  Differential  Element  for  Conservation  of  Momentum 
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Figure  2.4:  Control  Volume  for  Crack  Tip  Boundary  Condition 
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if  we   assume   a  unit   depth   in  the  z  direction.   Thus  the  momentum 
entering  the  left  face  per  unit  time  is: 

2 
(puw)u  =  pu  w  (2.14) 

The  momentum  leaving  the  right  face  per  unit  time  is: 

2 
2    d( pu   w)  ,  .  0  , , . 

pu  w  +  ■*-* ' dx  (2.15) 

r        dx 

The  momentum  leaving  the  top  and  bottom  faces  per  unit  time  is: 

2puv  dx  (2.16) 

2 
NOTE:    In  the  computer  coding,  Eq  2.16  was  incorrectly  coded  as  2pu  . 

For  the  conditions  considered  in  this  thesis,  it  was  determined 

that  this  error  was  inconsequential. 

By  defining  a  fanning  type  friction  factor  to  be: 

f  =  -? — 2  (2.17) 

. 5pu  v     ' 

where   r   is  the  shear  stress  on  the  top  and  bottom  sides,  we  can  find 

the  force  applied  as  a  result  of  this  shear  stress: 

Fg=  t   x  (area)  =  (|  pu2f)2dx  -  fpu2dx  (2.18) 

The  force  resulting  from  the  pressure  on  the  left  face  is: 

Pw  (2.19) 

The  force  resulting  from  the  pressure  on  the  right  face  is: 

Pw  +  aP*>dx  (2.20) 

dx 

Therefore,  the  sum  of  all  forces  being  considered  is: 

SF  =  (  -fpu2  -  3^)  )dx  (2.21) 

^       dx 

Then,   by   doing  a  force  and  momentum  balance  on  the  element,  we  find 

the  conservation  of  momentum  equation: 
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1  {  a^}  +  fpu2+  dif^K   2puv         (2.22) 
w    3x      r  dx 

dw  .    d(pu) 


2 . 2    Crack  Tip  Boundary  Condition 

To   derive   the   equations  for  the  crack  tip  a  slightly  different 
approach   is   used.    The   control  volume  used  is  shown  in  Figure  2.4. 

The  volume  is  assumed  to  be  a  wedge  with  V  =  t  ws  and  V  =  t  (ws  +  sw) , 

where   s   is   the   length   of   the   volume   and  w  is  the  crack  opening 
diplacement.   The  momentum  entering  the  control  volume  is: 

pu2w  (2.23) 

The  mass  leaving  the  volume  through  seepage  is: 

2pvs  (2.24) 

A  mass  balance  on  the  control  volume  yields: 

V^  =  pu2w  -  2pvs  -  pW  (2.25) 

For  the  conservation  of  energy  equation  we  first  need  to  find  the  flow 
work  and  energy  crossing  the  boundary.   The  flow  work  in  is: 

m  P  (2.26) 

P 

The  work  done  to  expand  the  volume  is: 

P^  (2.27) 

The  conservation  of  energy  expression,  therefore,  is: 


V  a^)=  (E  +  P)uw  +  fpu2sw  -(P  +  E)^7 
ot  at 


2sq"  -  2psv(E  +  P)         (2.28) 
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For   the   conservation  of  momentum  expression  we   need  to  know  the 
momentum  entering  and  leaving  the  volume  per  unit  time: 

(puw)u  -  (pvw)v  (2.29) 

and  the  net  force  acting  on  the  volume: 

(Pin-Pn)f  -  fpu2s  (2.30) 

where  P   is  the  pressure  found  at  the  crack  tip  node, 
n        r  r 

The   conservation  of  momentum  equation   for  the  crack  tip  volume  is 

therefore : 

v  aV     =  2  (p  ■  pn)  "  fpu  w  '  pu    at       (2-31) 


2 . 3  Crack  Mouth  Boundary  Condition 

At  the  mouth  of  the  crack,  the  high  pressure  gas  in  the  borehole 
expands  and  enters  the  crack.  At  any  time,  the  current  conditions  in 
the  borehole  and  in  the  crack  mouth  are  known.  The  task  is  to  predict 
the  flow  rate  into  the  crack  during  the  next  time  step.  The 
difficulty  is  that  that  the  pressure,  velocity,  and  density  of  the  gas 
entering  the  crack  are  not  the  same  as  the  conditions  in  the  crack  and 
that  this  state  cannot  be  chosen  arbitrarily. 

To  solve  this  problem,  we  pose  the  crack  mouth  conditions  as  a 
Riemann  problem.  In  the  Reimann  problem,  a  diaphragm  separates  two 
gases  at  different  states.  When  the  diaphragm  is  ruptured,  a  shock 
wave  and  rarefaction  wave  propagate  through  the  gases  (See  Chapter  A 
for  an  example  solution) .  The  solution  then  provides  a  complete 
definition  of  the  gas  conditions  passing  the  ruptured  diaphragm  (crack 
mouth) . 
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At  each  time  step  in  the  analysis,  the  conditions  in  the  borehole 

and   in   the  crack  mouth  are  known.   A  Riemann  problem  is  solved  using 

these  conditions.   The  solution  provides  the  state  of  the  gas  entering 

the  borehole   during  the  next  time  step.   This  is  then  repeated  for 
each  time  step. 

2  .4  Loss  Mechanisms 

Seepage  loss,  heat  transfer,  and  friction  expressions  follow  the 
approach  given  by  Nilson  et  al . ,  (1985).  The  lateral  seepage  loss 
into  the  wall  rock  is  estimated  with  a  one  dimensional  (normal  to  the 
rock  wall)  Darcy  flow-analysis: 

*-ii7rs^¥2  (232) 

where   P   is  the  local  gas  pressure  in  the  fracture,  P   is  the  ambient 

pore  pressure,  t  is  the  time  the  crack  face  has  been  exposed  to  gas,  k 
is  the  permeability,  $  is  the  porosity,  and  /i  is  the  viscosity. 

Heat  transfer  to  the  permeable  wall  rock  is  approximated  by  the 
expression: 

T  -  T 

conv     cond 
where   T   is   the   local   gas   temperature   in  the  fracture,  T   is  the 

ambient   temperature   of   the   rock,   R      is  the  resistance  to  heat 
r  conv 

transfer  by  convection,  and  R    ,  is  the  resistance  to  heat  transfer 
J  cond 

by  conduction.   R     is  estimated  from: 
J  conv 

4Pr2/  3 
R     =  T^?  (2.34) 

conv    fpC  u  v     ' 

P 
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where   Pr   is  the  Prantdl  number,  C   is  the  constant  pressure  specific 

P 

heat  of  the  gas,  and  u  is  the  local  velocity  of  the  gas. 

The   resistance   to  heat   transfer  through  conduction  based  on  a 
one -dimensional  lateral  diffusion  analysis  is: 

R      .  (Q  I        t  )  (2  35) 

cond     k 


wh 


ere   a   is   the  thermal  diffusivity,  k   is  the  thermal  conductivity, 

and  t  is  the  time  the  wall  is  exposed  to  the  gas. 

The   frictional   effects   for  both  turbulent  and  laminar  flow  are 
incorporated  into  a  single  expression: 

f  =  ^  +  0.055  (  -  )0-472  (2.36) 

Re  w 

where  Re  is  the  Reynolds  number  (  ^—   )  ,  /i  is  the  gas  viscosity,  and  e 

is  the  fracture  roughness.  During  low  Re  for  laminar  flow,  the  first 
term  is  predominant.  During  high  Re  for  turbulent  flow  the  second 
term  becomes  predominant  with  a  smooth  transition  between  the  two 
extremes . 

A  Sutherland  viscosity  model  used  for  air  is  used  to  estimate  the 
viscosity  of  the  gas  (Schlichting,  1968): 

U.  _   rl  ,3/2  V  Sl  3 

u    lT  J     T  +  S,  K*>*'J 

^o     o  1 


wh 


ere  \x     denotes  the  viscosity  at  the  reference  temperature  T  ,  and  S, 


is  a  constant  which  for  air  has  the  value  S,=  110  K, 
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2 . 5  Finite  Difference  Approximation 

Since   the   equations   of  gas  dynamics  are  based  on  conservation 

laws ,  researchers  in  fluid  mechanics  have  often  found  it  useful  to  use 

a   form  of  the   equations,  called  the  divergence  form,  which  clearly 

displays  the  conserved  quantities  of  mass,  momentum,  and  energy  (Ames, 

1977).    It  would   therefore  seem  reasonable  to  try  to  preserve  these 

conservation  properties   in  the  finite  difference  approximations.   A 

system  of  equations: 

aw.  +  3JLIw)=0  3 

dt  3x 

where   w  is  a  vector  funtion  (of  x  and  t  with  n  components)  and  f  is  a 

nonlinear  vector   function   (of  the  vector  w  with  n  components)  is 

called  a  system  of  conservation- law  form.   The  original  suggestion  of 

Lax  was   a   staggered  scheme  which,   when  applied  to  Equation  2.38, 

gives : 

k  '  "i,J+r  I  ("i+i,j  +"i-i,j>  > 

2h  <  fi+i,r£i-i,j  '  -°  <2-39> 

where   i   is   the  node  position,  j  is  the  current  time,  k  is  the  time 
step,  and  2h  is  the  distance  between  the  i-1  and  i+1  nodes. 
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Writing  conservation  Equations  2.5,  2.11,  and  2.21  in  this  form  we  get 
(all  terms  on  the  right  hand  side  of  the  following  equations  are  at 
the  current  , or  j ,  time): 

i  ,  ,  L  At  ,  (puw)i-r  (^uw)i-n 

pi,j+l  =  2  (pi+l+  pi-l}  +  w  {        Ax 


-  (2v  +  ^  )p.  }         (2.40) 


pui,j+l  "  2  (pUi+l  +  pUi-l} 


+  A£  (P  +  (piQu).^-  (P  +  (pu)u).+1 
w  Ax 

-fpu  u  -  2puv  -    pu     —    }  (2.41) 

ot 


Ei,j+1=2  <Ei+l+Ei-l> 


^t    (uw(E  +  P))i.1-  (uw(E  +  P))i+1 
w  Ax 

•(P  +  E)^7  -  2v(E  +  P)  -  2q"  }  (2.42) 

ot 


2.5.1  Stability 

The  nodal  staggering  enables  central  space  difference  and  forward 
time   difference   to  be  used  without  developing  instability.   These 

staggered  schemes   are  usually  stable  if  (  a  +  v   )  <  —  where  a  is 

the  local  sound  speed  and  v  is  the  local  gas  velocity. 

In  essence,   the   replacement  of  w.  .  by  the  average  of  neighboring 

•*•  >  J 

values  in  Equation  2.39  has  the  effect  of  introducing  a  dissipative 
term.  The  Lax  scheme  is  thus  said  to  have  numerical  dissipation,  or 
numerical  viscosity  (Press  et  al . ,  1986). 
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For  an  intuitive  explanation  for  stability,  refer  to  Figure  2.5. 
The  quantity  u.  .  ,  is  computed  from  information  at  points  i-1  and  i+1 

J.  ,  JTl 

at   time  j  (or  t) .   In  other  words,  x.  , and  x.  .are  the  boundaries  of 

the   spatial   region  that   is   allowed  to  communicate  information  to 

u.   ,,.    Recall   that   in  a  continuum  wave   equation,   information 
l,j+l  H 

actually  propagates  with  a  maximum  velocity  v.   If  the  point  u.  .  ,  is 

outside  of  the  shaded  region  in  Figure  2.5,  then  it  requires 
information  from  points  more  distant  than  the  differencing  scheme  will 
allow.  The  lack  of  that  information  will  result  in  the  instability. 
Therefore,  At  cannot  be  made  too  large. 

The  goal  of  recent  work  on  finite  differencing  schemes  is  to  make 
the  discontinuities  from  shocks  and  contact  surfaces  more  accurate. 
As  will  be  discussed  later,  when  the  energy,  mass,  and  momentum  losses 
to  be  considered  are  added  to  the  conservation  equations,  these  sharp 
discontinuities  will  not  be  present.  Because  of  this,  the  Lax  method 
should  give  sufficiently  accurate  results. 

2.6  Burn  Model 

The  equations  used  for  this  procedure  are  presented  from  Schatz 
and  Hanson  (1986)  and  Mniszewski  and  Napadensky  (1985).  The  mass  of 
the  propellant  burned,  M(t) ,  as  a  function  of  time  may  be  expressed 
as : 

M(t)  =  M   f(t)  [  1  +  K  -  K  f(t)  ]        (2.43) 

where  M   is  the  initial  mass  of  propellant  present  in  the  borehole,  K 

is   the   form  coefficient,  and  f(t)  is  the  fraction  of  web  burned  as  a 
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Figure  2.5:  Graphical  Representation  of  Stability 
(After  Press  et  al . ,  1986) 
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function  of  time.  The  fraction  of  web  burned,  f (t) ,  is  expressed  by 
the  differential  equation: 

d    f(t)  _  2  B  P(t)  n  n    ... 

at  D  {*.*<*) 

where  B  is  the  burn  rate  constant,  P(t)  is  the  current  borehole 
pressure,  D  is  the  web  thickness,  and  n  is  the  burn  rate  exponent  (See 
Figure  2.6).  The  mass  of  gas  in  the  borehole  as  a  function  of  time 
is : 

"bh"11^  -  Mck(t>  (2-45) 

where   M(t)   is   the   mass  of  gas  in  the  borehole  if  there  was  no  loss 

into   the   crack,   and  M  ,  (t)  is  the  mass  lost  to  the  crack.   The  net 

ck 

mass   leaving  the  borehole  is  accumulated  in  the  subroutine  GAS_SOLVE. 

The  borehole  volume  available  for  gas  expansion  is  given  by: 

V,,=  V.  +  M  [  ~   +  f(t)  '  l    }  (2.46) 

bh    l    oL  p  p  J 

*p     ^s 

where  V.   is   the   initial  volume   of  air  surrounding  the  propellant 

canister,  p        is   the  pack  density,  and  p      is  the  solid  density.   The 
P  s 

change  of  volume  due  to  expansion  of  the  borehole  was  ignored. 
The  pressure  in  the  borehole  is  then  found  using: 


P(t)  =        Mbh(t)  F  (2.47) 

Vbh^)  *  t  Mbh(t)  C  1 

where  F  is  the  propellant  energy  density,  and  C  is  the  covolume. 
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a.   Unburned  Grain 


I —    Slivers     — i 


b.   Burning  Grain 


Figure  2.6:  Web  Thickness  and  Route  of  Burning  Through  a  Mult-Perf 
Progressively  Burning  Grain  (After  Mniszewski  and 
Napadensky,  1985) 
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CHAPTER  3 
GAS  SOLVE 

This  chapter  will  describe  the  features  available  to  the  user  to 
set  up  a  gas  flow  problem.  The  Gas_Solve  subroutine  was  incorporated 
into  CRACKER,  a  finite  element  code  used  for  modeling  dynamic 
propagation  of  a  discrete  crack.  CRACKER  uses  the  finite  element 
method  to  model  the  continuum  and  automatically  remeshes  as  the  crack 
propagates.  For  information  on  CRACKER  and  its  uses,  see  Swenson 
(1985).  Also,  a  description  of  the  solving  routine,  Gas_Solve,  will 
be  given. 

3 . 1  Gas  Parameters 

Most  of  the  information  needed  to  set  up  a  gas  flow  problem  is 
found  in  the  gas  parameters  page.  A  sample  (used  in  Case  5,  Chapter 
5),  is  shown  in  Figure  3.1.  Every  value,  except  the  current  time 
step,  can  be  interactively  changed  by  the  user. 

The  user  has  the  option  to  choose  a  compressible  gas  model,  or 
to  specify  a  known  pressure  distribution  along  the  crack  face.  If 
the  gas  model  is  chosen,  the  user  is  prompted  to  choose  the  solution 
method  to  be  used.  The  Lax  method  is  currently  the  only  solution 
method  in  use,  and  is  illustrated. 

The  initial  gas  conditions  in  the  crack  are  specified  to  be 
atmospheric  conditions.  The  gas  constants  are  chosen  to  be  that  of 
nitrogen.  If  the  burn  model  option  is  chosen,  the  values  describing 
the  borehole  and  charge  volumes,  and  propellant  properties  must  be 
specified.   Another  option  is  to  specify  pressure  and  temperature  time 
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Pressure  in  cracks =   GAS  MODEL 

Method  used  to  solve  gas  tlo  w  equations  =    LAX 

Initial  crack     pressure <  t  =  0  ) =  0.10000E+06 

Initial  crack  t  e  m  P  e  r  a  t  u  r  e  <  t  =  9  > =  0.30000E+03 

Specific  heat  ratio  <k> =  0.14000E+01 

Ideal  gas  constant  (R) =  0.£3688E+83 

To  define  borehole  conditions  we  will  use  A  BURN  MODEL 

Borehole  u  o  1  u  m  e =  8.78548E-02 

Solid  density =  0.16500E+04 

Pack  density  =  y  ^6w9SE,t'33 

Form  coefficient. =  8.62888E  +  98 

Burn  rate  coefficient =  0.539SSE-O7 

Burn  rate  exponent =  9  8800  0E+00 

Covolume =  9.16099E-0E 

W  e  b  thickness =  8.76288E-83 

Pre.  pell  ant  energy  density =  0.10337E+07 

Charge  volume =  0  19635E-9d 

Use  friction  model V  E  S 

Friction  factor  <  F  > =  8.18008E-01 

Surface  roughness =  0  29999E-0  4 

Use  seepage  model Y  E  S 

Permeability =  0  10808E-14 

Porosity =  0.E0000E-01 

Reference  viscosity =  8.18000E-04 

Far  field  Pressure(rock) =  0.18888E+8S 

Sutherland  constant =  0.11800E  +  03 

Reference  temperature  for  SI =  0.38080E  +  03 

Use  heat  transfer  model YES 

Far  field  temperature(rock) =  0.39999E+93 

Conductivity =  0  S8008E  +  81 

Rock  thermal  diffusiuity =  0.40008E-05 

Time  integration  factor =  0.98800E  +  80 

Damping  factor =  0.80000E  +  00 

Delta  x =  0.30000E-01 

Default  crack     width =  0  10000E-82 

Current  (gas)  time  step =  8.1546SE-04 


Figure  3.1:  Sample  of  Gas  Parameters  Page 
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histories  for  the  borehole.  Examples  of  these  are  shown  in  Figure  5.3 
(Case  2,  Chapter  5).  These  plots  are  input  interactively  in  the 
boundary  conditions  page  found  in  CRACKER. 

If  the  friction  model  is  used,  then  Equation  2.36  will  be  used. 
If  not  used,  then  the  friction  factor  must  be  specified.  If  the 
factor  is  zero,  then  no  friction  effects  will  be  seen.  If  the  factor 
is  nonzero,  then  the  specified  constant  friction  factor  will  be  used. 

If  the  seepage  model  is  chosen,  then  the  values  describing  the 
rock  and  far  field  conditions  must  be  specified.  If  this  option  is 
not  taken,  no  seepage  effects  will  be  seen. 

If  the  heat  transfer  model  is  chosen,  then  the  values  describing 
the  rock  and  far  field  conditions  must  be  specified.  If  this  option 
is  not  taken,  no  heat  transfer  effects  will  be  seen. 

The  time  integration  factor  is  multiplied  by  the  calculated  time 
step  to  result  in  the  time  step  actually  used.  This  reduction  factor 
should  be  kept  as  close  to  one  as  possible  to  prevent  smearing  of  the 
solution. 

The  damping  factor  should  also  be  kept  as  close  to  a  value  of 
one  as  possible.  It  was  found  that  in  the  early  stages  of  a  problem, 
the  model  became  unstable  when  friction  was  introduced.  To  promote 
stability,  the  damping  was  introduced.  For  the  interior  cells,  the 
damping  results  in  P  %  of  the  new  value  to  be  added  to  (1-P)  %  of  the 
old  value.  This  sum  will  be  the  new  value.  At  the  inlet  cell  (P-.3)  % 
is  used.  For  this  example,  damping  was  0.8.  As  will  be  seen  in  the 
verifications  chapter,  the  damping  has  a  minimal  effect. 

The  delta  x  (node  spacing)  was  chosen  to  be  .03  meters. 
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If,  during  the  course  of  the  analysis,  the  crack  opening 
displacement  was  found  to  be  less  than  the  default  crack  width,  the 
diplacement  was  set  to  the  default.  It  was  found  that  when  friction 
was  once  again  introduced,  instabilities  were  a  problem  at  the  initial 
stages  of  a  problem  when  the  crack  was  short  and  closed. 

3.2  GAS_SOLVE 

The  flow  path  and  logic  used  in  GAS_SOLVE  is  shown  in  Figure 
3.2.  The  GAS_SOLVE  routine  is  inside  a  larger  loop  in  CRACKER.  If 
new  crack  face  surface  tractions  are  needed,  GAS_SOLVE  is  entered. 
The  first  time  GAS_SOLVE  is  used,  arrays  are  initialized  which 
includes  initializing  gas  conditions  in  the  crack.  Then  a  loop  over 
all  nodes  is  entered.  Inside  this  loop  all  mass,  momentum,  and  energy 
losses  are  taken  care  of  on  a  node  by  node  basis.  Once  all  necessary 
calculations  are  finished,  a  comparison  of  the  current  time  in 
GAS_SOLVE  (GASJICURR)  and  CRACKER  (TCURR)  is  performed.  If  GAS_TCURR 
is  greater  than  TCURR,  then  an  interpolation  through  time  is 
performed,  and  GAS_SOLVE  is  exited.  Otherwise,  GASJICURR  is 
incremented  by  a  time  step  and  the  process  is  repeated. 


27 


IF   t=0.  Initialize 


Calculate   crack   opening   dlspl.   &   vel. 


Save  old  data 


Calculate   time   step 


loop  on  alt   cells 


Calculate  losses 


Calculate   gas   conditions 


Update   gas/crack   face   exposure   tine 


Check   if   finished 


ND 


YES 


Match  time  found  in  CRACKER 


Return  to  CRACKER 


t+dt 


Figure  3.2:  The  GAS_SOLVE  subroutine 
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CHAPTER  4 
VERIFICATION  PROBLEMS 

4 . 1  Burn  Model 

Closed  bomb  tests  were  conducted  by  the  Gas  Research  Institute 
(Mniszewski  and  Napadensky,  1985)  to  determine  the  properties  of 
various  propellants,  particularly  the  burn  rate  constant  and  burn  rate 
exponent.  The  original  closed  bomb  test  results  are  shown  in  Figure 
4.1.  Using  these  results,  the  burn  model  used  in  this  thesis  was 
verified.  The  experimentally  determined  values  for  the  burn  rate 
exponent  and  burn  rate  constant,  in  addition  to  the  other  quantities 
shown  in  Table  4.1,  were  used  in  the  burn  model.  The  verification 
results  are  shown  in  Figure  4.2.  Although  there  is  a  14  percent  error 
when  comparing  the  experimentally  obtained  peak  value  to  the 
calculated  peak  pressure,  the  time  at  which  the  peak  pressure  occurs 
is  in  close  agreement.  All  of  the  differences  involved  are  probably 
the  result  of  errors  encountered  when  estimating  the  constant  values 
from  the  data. 

4.2  Shock  Tube  Problem 

The  differential  equations  were  cast  in  finite  difference  form 
by  using  the  Lax  method.  It  was  thought  that  with  the  simplicity 
gained  by  using  this  method,  the  first  order  accuracy  the  Lax  method 
delivers  could  be  accepted.  To  get  a  feeling  for  the  accuracy  lost,  a 
typical  shock  tube  test  problem  was  run. 

The  shock  tube  is  a  useful  research  tool  in  which  normal  shock 
waves  are  generated.   They  are  used  in  the  study  of  materials  that  are 
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Figure  4.2:  Numerical  Results  of  Closed  Bomb  Test 
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Propellant 

Amount  of  propellant 

Web  thickness 

Volume 

Energy  density 

Covolume 

Burn  rate  constant 

Burn  rate  exponent 

Solid  density 

Pack  density 


NQ/M 

10  grams 

0.001118*1.15  m 

8.194E-5  m3 

1052. 5E3  J/Kg 

0.001085  m3/kg 

4.0653E-7  m/sec/Pa 

0.6495 

1665  kg/m3 

750  kg/m3 


Table  4.1:  Burn  model  verification  parameters 
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subjected  to  the  extreme  pressure  and  temperature  conditions  found 
behind  a  shock  wave  (John,  1984) .  This  is  accomplished  by  separating 
a  high  pressure  gas  from  a  low  pressure  gas  using  a  diaphragm. 
Typical  pressure,  density,  and  velocity  profiles  shortly  after  rupture 
of  the  diaphragm  are  shown  in  Figure  4.3  (Sod,  1978).  when  the 
diaphragm  is  ruptured,  a  discontinuity  in  gas  properties  exists,  which 
is  characteristic  of  a  normal  shock.  This  normal  shock  then  moves 
into  the  low  pressure  side,  with  a  series  of  expansion  waves 
propagating   into   the  high  pressure  gas.   Points  x,  and  x«  represent 

the  location  of  the  head  and  tail  of  the  expansion,  or  rarefaction, 
wave  moving  to   the   left.    The  point  x^   is  where   the   contact 

discontinuity  occurs.   The  pressure  and  velocity  are  continuous  across 

the   contact   surface.    However,   the  density  and  temperature  are  not 

continuous   across  a  contact  surface.   Point  x.  is  the  location  of  the 

4 

shock  wave   moving   to   the   right.    Across   the   shock,   all  of  the 

quantities  (p,T,P,U)  will  be  discontinuous. 

In  this  test  of  the  model,  the  following  conditions  were  used: 

p,=  1.0        P1=  1.0       u,=  0.0 

P5=  0.125      P5=  0.1       u5=  0.0 

The  ratio  of  specific  heats  was  1.4  and  Ax  =  .002  meters.  The 
analytic  solution  of  the  hyperbolic  system  of  conservation  laws 
(Riemann  problem)  was  found  using  a  Brent  iteration  scheme  (see 
Halter,  1985).  Figures  4.4,  4.5,  and  4.6  compare  the  pressure, 
density,  and  velocity  distributions  obtained  using  the  Lax  and  Brent 
methods . 
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Density 

Pressure 

Density 
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Figure  4.3:  Shock  Tube  t  >  0  (After  Sod,  1978) 
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Figure   4.4:    Pressure   profile   of  Riemann  problem  using  Lax  method 


35 


ELRPSED  TIME=   0.  1443E+00 


O      1.200CH- 
X 


1.0000 


0.  B000-- 


(—       0.  6000- 
M 

Z 
UJ 
Q      0.  4000- ■ 


0.  2000-- 


0.0000- 


-t- 


■+■ 


-V 


•+■ 


0.0000         0.2000         0.4000         0.6000         0.9000  1.0000 


DISTANCE 


(X10  ) 


Figure  4.5:  Density  profile  of  Riemann  problem  using  Lax  method 
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Figure  4.6:  Velocity  profile  of  Riemann  problem  using  Lax  method 
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There  is  relatively  close  agreement  across  the  rarefaction. 
However,  the  corners  at  the  endpoints  of  the  rarefaction  are  rounded. 
There  is  not  as  close  agreement  at  the  contact  surface  and  shock.  The 
reader  should  keep  in  mind  that  the  purpose  of  this  model  is  to 
predict  gas  flow  under  conditions  that  are  definitely  not  isentropic. 
Isentropic  conditions  are  assumed  for  shock  tubes.  Because  of  this, 
the  inability  of  the  Lax  method  to  accurately  predict  sharp 
discontinuities  should  not  effect  the  results  greatly. 

4. 3  Simplified  Crack  Problem 

Although  this  section  does  not  describe  a  verification  problem, 
it  does  give  some  insight  into  the  effects  of  some  of  the  variables 
involved.  The  results  to  be  presented,  in  addition  to  others,  were 
used  to  debug  the  coding,  and  also  to  begin  to  understand  gas  flow  in 
a  crack  before  large  problems  were  attempted.  Different  ideas  could 
be  tried  and  tested  with  these  problems  using  minutes  of  computer 
time.  Typical  large  problems  ran  in  excess  of  15  CPU  hours  on  a  H800 
Harris  mini  -  computer . 

A  100  meter  crack  with  uniform  cross-sectional  area  was  used. 
The  borehole  temperature  and  initial  crack  temperature  was  initially 
1000  K.  The  initial  borehole  pressure  was  40MPa  and  the  initial  crack 
pressure  was  1  MPa .  The  cases  considered  are  displayed  in  Table  4.2. 
The  variables  f,  v,  and  q  represent  friction,  seepage,  and  heat 
transfer. 

The  purpose  of  Cases  1  &  2  was  to  determine  the  effect  of  the 
time  step  reduction  factor  and  damping.  Figures  4.7  and  4.8  show  the 
pressure   and  velocity  results  for  Cases  1  &  2.   The  flow  in  the  crack 
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Initial 

Borehole  Crack  time 

Pressure  Pressure  step 

Case     (MPa)  (MPa)  reduction 

1  40E6      1E6  1 

2  40E6      1E6  0.9 

3  40E6      1E6  0.9 


damping   A  x 

factor     (m)  losses 

0.99    0.205  '  none 

0.8     1.01  none 

0.8     1.01  f,v,q 


Table  4.2:  Verification  Problems 
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Figure  4.7:  Profiles  for  Verification  Case  1 
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Figure  4.8:  Profiles  for  Verification  Case  2 
40 


is  supersonic,  and  the  rarefaction  wave  is  moving  into  the  crack 
mouth.  This  causes  the  reduction  from  the  borehole  pressure  at  the 
crack  entrance.  Although  the  shock  wave  isn't  as  distinct  in  Case  2, 
the  effect  of  the  pressure  wave  is  still  present.  Figure  4.9  shows 
the  pressure  and  velocity  profiles  for  Case  3.  The  velocity  remains 
subsonic  and  the  pressure  decreases  rapidly  from  the  crack  mouth.  The 
smoothing  of  the  shock  wave  with  the  introduction  of  friction,  heat 
transfer,  and  seepage  is  typical  of  what  is  expected  when  propagating 
a  crack. 
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CHAPTER  5 
APPLICATIONS 

In  this  chapter  we  address  the  question  of  possible  fracture 
length  for  realistic  burn  conditions.  To  do  this,  we  use  a  quarter 
model  of  a  borehole  with  symmetry  conditions  imposed.  Therefore,  the 
model  represents  a  borehole  with  four  radial  cracks .  Several 
different  borehole  pressures  and  loading  rates,  using  either  a  burn 
model  for  a  propellant  or  specified  pressure  time  histories,  are  used. 

The  model  geometry,  shown  in  Figure  5.1,  was  used  in  Cases  3  & 
4.  This  geometry  was  used  to  minimize  the  required  computer  run  time 
after  it  was  realized  long  fractures  were  not  going  to  be  seen.  The 
other  problems  used  a  square  geometry  that  extended  70  meters.  The 
material  properties  for  sandstone  are  shown  in  Table  5.1.  The  rock  and 
gas  values  needed  by  the  gas  model  are  shown  in  Table  5.2.  No  dynamic 
critical  stress  intensity  data  is  known  for  sandstone,  so  the 
estimated  dynamic  curve  was  based  on  the  static  value.  The  curve  used 
is  shown  in  Figure  5.2. 

A  borehole  diameter  of  0.2  meters  was  assumed  with  an  initial 
0.146  meter  fracture  extending  at  45  degrees  to  the  X  and  Y  axis. 
Since  symmetry  was  present  in  these  cases,  this  actually  represents 
four  fractures  in  a  radial  pattern.  The  outer  boundary  was  at  least 
35  meters  from  the  borehole  in  all  cases  to  prevent  interference  from 
a  reflected  stress  wave. 

To  set  up  the  problem,  the  initial  in-situ  stresses  were  first 
applied  uniformly   to   the  mesh.   For  all  cases,  the  in-situ  stresses 
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BOREHOLE 


Figure  5.1:  Model  Geometry 


Figure  5.2:  Assumed  Dynamic  Stress  Intensity  for  Sandstone 
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E  =  24.82E9  Pa 
v   =  0.22 
p   =  2600  Kg/m3 
KIc  =  1.978E6  PaTm 


c.=  3048  m/sec 


c  =  1806  m/sec 


c  =  1651  m/sec 
R 


Table  5.1:  Material  properties  for  sandstone 
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Initial  crack  pressure  100  KPa 

Initial  crack  temperature  300  deg  K 

Specific  heat  ratio  1.4 

Ideal  gas  constant  296.8  J/Kg/K 
Borehole  volume (qtr  symmetry)   0.007854  m3 

Charge  volume(qtr  symmetry)  0.0019635  m3 

Solid  density  1650  J/Kg 

Pack  density  700  J/kg 

Form  coefficient  0.62 

- 1     8 

Burn  rate  coefficient  0.53968E-7  m  sec   Pa 

Burn  rate  exponent  0.8 

Covolume  0.001  m3/Kg 

Web  thickness  0.000762  m 

Propellant  energy  density  1038. 7E3  J/Kg 

Surface  roughness  5-20/im 

Permeability  0.1E-14  m2 

Porosity  0.02 

Far  field  pressure  100  KPa 

Far  field  temperature  300  deg  K 

Conductivity  2  J/m/K/sec 

Thermal  diffusivity  0.4E-5  m2/sec 


Table  5.2:  Parameters  used  in  gas  model 
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were  a   =  a   =   -17.23MPa.   No  overburden  stress  was  assumed.   A  short 
x    y 

test  using  a   =   -6.89MPa  found  that  there  was  a  10-14%  difference  in 

crack  length  and  KT  when  o     was  ignored.   Before  beginning  the  dynamic 

analysis,  static  equilibrium  was  achieved  by  introducing  a  small 
amount  of  damping  and  allowing  the  problem  to  run  until  there  was 
negligible  change  in  node  displacement.  This  allows  the  stresses 
around  the  borehole  to  come  to  equilibrium  before  proceeding  with  the 
pressure  loading.  For  the  cases  presented  here,  equilibrium  occured 
at  t  =  0.002  sec.  The  reader  should  keep  this  in  mind  when  evaluating 
time  plots  since  the  actual  dynamic  problem  starts  at  2  msec. 

5 . 1  Case  1  -  Specified  Borehole  Pressure  With  No  Losses  In  The 
Fracture 

Case  1  will  be  used  as  the  reference  case  with  which  other 
results  will  be  compared.  Although  a  no-loss  case  is  not  realistic, 
it  does  give  the  easiest  case  for  comparison.  The  borehole  pressure 
was  ramped  to  103.4  MPa  and  the  borehole  temperature  was  ramped  to 
3000  deg.  K.  These  time  histories  are  shown  in  Figure  5.3.  Figure  5.4 
shows  the  displaced  mesh  results  at  12  msec.  At  this  time  the  crack 
is  about  8.5  meters  long  and  the  crack  opening  displacement  at  the 
borehole  is  approximately  0.01  meters.  Plots  of  crack  opening  and 
crack  length  are  shown  in  Figures  5.5  and  5.6.  The  crack  tip 
velocity,  Figure  5.7,  reaches  the  limiting  crack  velocity  of  0.5  of 
the  Rayleigh  wave  speed  during  propagation.  This  limit  is  imposed  by 
the  numerical  solution  scheme.  In  reality,  bifurcation  would  probably 
occur   at   this   large  loading.   The  significant  point  is  that  without 
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Figure  5.3:  Borehole  Time  Histories 
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Figure  5.4:  Displaced  Mesh  Results  for  Case  1  at  12  msec 
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losses,  the  limit  on  crack  velocity  is  due  to  the  rock  response.  The 
gas   is   attempting  to  over-drive   the  crack.   Plots  of  K,„  and  KTT_. 

during  propagation  are  shown  in  Figures  5.8  and  5.9.  Comparing  the 
crack  tip  velocity  to  the  gas  velocity  shown  in  Figure  5.10,  it  can  be 
seen  that  the  pressure  wave  has  reached  the  crack  tip  and  is  partially 

reflected.    This  will   result  in  the  large  KT_  and  the  fact  that  the 

crack  is  driven  hard  enough  to  probably  bifurcate.  Temperature, 
pressure,  density,  and  mach  number  plots  are  shown  in  Figures  5.11 
through  5.14.  Note  that  the  gas  flow  is  supersonic  for  a  significant 
portion  of  the  crack  length.  Since  the  pressure  wave  has  reached  the 
crack  tip,  all  the  property  values  are  relatively  constant  along  the 
crack  face.  It  will  be  seen  that  when  losses  are  introduced,  there 
are  large  variations  between  property  values  at  the  crack  mouth  and 
crack  tip. 

To  illustrate  that  the  gas  pressure  is  being  correctly  applied 
to  the  crack  face,  Figure  5.15  shows  the  normal  stress  along  a  line 
parallel  to  the  crack.  The  normal  stress  distribution  is  similar  to 
the  gas  pressure  in  the  fracture. 

An  interesting  result  of  this  analysis  is  shown  in  Figure  5.16. 
From  this  plot,  the  net  energy  that  has  left  the  borehole  up  to  a 
given  time   can  be   found.    In  future  reference,  this  plot  will  be 

called  the  crack  entrance  energy  plot.  At  6  msec,  about  0.96x10  J 
entered  the  crack  and  has  driven  the  fracture  about  3.0  meters.  It 
should  be  kept  in  mind  that  some  of  this  energy  has  not  been  used  to 
drive  the  crack,  but  has  instead  been  used  to  increase  the  kinetic  and 
internal   energy  of  the  gas  in  the  fracture.   This  stored  energy  will 
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Figure  5.7:  Crack  Velocity  for  Case  1 
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Figure  5.8:  KID  During  Propagation  for  Case  1 
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Figure  5.9:  KIID  During  Propagation  for  Case  2 
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Figure   5.10:    Gas  Velocity   for   Case   1   at   6.5,    9.5,    and   12   msec 
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Figure   5.11:    Gas  Temperature   for  Case   1   at   6.5,    9.5,    and  12   msec 
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Figure  5.12:  Gas  Pressure  for  Case  1  at  6.5,  9.5,  and  12 


msec 
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Figure  5.13:  Gas  Density  for  Case  1  at  6.5,  9.5,  and  12  msec 


Figure  5.14:  Mach  Number  for  Case  1  at  6.5,  9.5,  and  12  msec 
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Figure  5.15:  Normal  Stress  along  Crack  Face  at  12  msec 
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Figure  5.16:  Crack  Entrance  Energy  for  Case  1 
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later  become  available  to  further  drive  the  crack.  Figure  5.17  shows 
the  work  done  to  move  the  boundaries  of  the  problem.  Although  the 
work  done  to  expand  the  borehole  is  included,  its  contribution  should 
be   minimal.    At  t=0 ,  this  value  is  initialized  by  using  the  combined 

initial  values   of   strain  and  kinetic   energy   (E  n»  45.7x10  J). 

Subtracting   this  initial  value  from  the  value  shown  in  Figure  5.17  at 

6  msec  will  give  us  the  boundary  work  done  of  0.3x10  J.  If  we 
subtract   the   energy  doing  boundary  work  from  the  supplied  energy,  an 

approximation  of  the  energy  present  in  the  gas  of  0.93x10  J  can  be 
found.  Comparing  this  value  to  the  calculated  value  of  gas  energy  at 
6  msec,  shown  in  Figure  5.18,  we  can  find  there  is  a  3  per  cent 
difference.  The  majority  of  this  error  is  caused  by  the  inclusion  of 
the  borehole  expansion  work.  Even  with  this  error,  the  conservation 
of  energy  is  confirmed  to  hold  true  in  this  analysis. 


5.2  Case  2    -  Specified  Borehole  Pressure  With  Friction.  Heat  Transfer . 
And  Seepage  Included 

For  case  2,  the  borehole  pressure  and  temperature  were  again 
ramped  to  the  maximum  values,  as  in  case  1.  At  12  msec,  the  crack 
mouth  opening  and  crack  length  are  about  0.85  cm  and  4  meters, 
respectively.  Plots  of  crack  opening  and  crack  length  are  shown  in 
Figures  5.19  and  5.20.  In  contrast  to  case  1,  the  crack  speed  is 
controlled  by  the  gas  flow  into  the  crack.  The  crack  is  not  being 
overdriven.  The  displaced  mesh  is  shown  in  Figure  5.21.  From  the 
crack  entrance  energy  plot  shown  in  Figure  5.22,  we  can  see  that  about 
67%   of   the   energy  required  in  Case  1  is  used  to  run  a  crack  that  is 
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Figure  5.17:  Energy  Required  to  Perform  Boundary  Work  for  Case  1 

(Initialized  at  45.7E7  J) 
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Figure  5.18:  Energy  Present  in  Gas  for  Case  1 
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Figure   5.19:    Crack  Opening  Displacement   for  Case   2   at   12  msec 
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Figure  5.20:  Crack  Length  for  Case  2 
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Figure  5.21:  Displaced  Mesh  Results  for  Case  2  at  12 


msec 
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Figure  5.22:  Crack  Entrance  Energy  for  Case  2 
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about  47%  the  length  found  in  Case  1.  This  difference  in  fracture 
length  is  caused  by  heat  transfer,  seepage,  and  friction  losses. 
Figures  5.23  and  5.24  show  heat  transfer  and  seepage  losses.  Note  the 
decreasing  values  as  a  function  of  time.  The  loss  mechanisms 
introduced  in  this  case  tend  to  convert  more  of  the  kinetic  energy 
present  at  the  crack  mouth  to  internal  energy.  This  results  in  the 
large  variation  of  gas  property  values  from  the  crack  mouth  to  the 
crack  tip.  This  can  be  seen  in  Figures  5.25  through  5.29  where  gas 
properties  are  shown.  While  the  pressures  and  temperatures  at  the 
crack  mouth  are  higher  than  those  found  in  Case  1,  the  velocities  are 
lower  illustrating  this  energy  conversion.  Figure  5.30  shows  a  plot 
of  Kjn  during  propagation. 

5 . 3  Case  3  -  Specified  Borehole  Pressure  With  Friction  In  The  Fracture 
Included  (  e  =  20/xm  ) 

This  case  is  identical  to  the  previous  case  except  that  the 
effects  of  seepage  and  heat  transfer  are  not  included.  By  comparing 
these  results  to  Case  2,  we  can  determine  the  individual  effect  of 
friction,  and  the  combined  effect  of  heat  transfer  and  seepage. 

At  12  msec  the  crack  length  is  about  7  meters.  Plots  of  crack 
length  and  KTn  are  shown  in  Figures  5.31  and  5.32.   Again,  the  limit 

on  crack  velocity  is  a  result  of  reduced  gas  flow  into  the  crack.  The 
displaced  mesh  is  shown  in  Figure  5.33.  Figure  5.34  shows  the  crack 
tip  velocity  during  propagation.  Figure  5.35  shows  a  crack  entrance 
energy  plot  for  Case  3.  The  slopes  of  this  graph  and  the  crack 
entrance  energy  plot  for  Case  2  are  similar.  This  would  imply  that 
about   the   same  percentage  of  the  energy  available  in  Case  3  is  being 
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Figure  5.23:  Heat  Transfer  Results  for  Case  2  at  6.5,  9.5, 

and  12  msec 


Figure  5.24:  Seepage  Results  for  Case  2  at  6.5,  9.5,  and  12  msec 
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Figure   5.25:    Gas   Pressure   for  Case   2   at   6.5,    9.5,    and   12  msec 
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Figure  5.26:  Gas  Velocity  for  Case  2  at  6.5,  9.5,  and  12  msec 
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Figure  5.27:  Gas  Temperature  for  Case  2  at  6.5,  9.5,  and  12  msec 
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Figure  5.28:  Gas  Density  for  Case  2  at  6.5,  9.5,  and  12  msec 
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Figure  5.29:  Mach  Number  for  Case  2  at  6.5,  9.5,  and  12  msec 
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Figure  5.30:  KID  During  Propagation  for  Case  2 
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Figure  5.31:  Crack  Length  for  Case  3 
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Figure  5.32:  KID  During  Propagation  for  Case  3 
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Figure   5.33:    Displaced  Mesh  Results   for  Case   3   at   12 


msec 
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Figure  5.34:  Crack  Velocity  for  Case  3 
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Figure  5.35:  Crack  Entrance  Energy  for  Case  3 
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converted  to  kinetic  and  internal  energy.  This  should  result  in 
similar  gas  property  plots  as  can  be  seen  in  Figures  5.36  through 
5.40. 

When  comparing  Cases  2  &  3,  it  can  be  seen  that  seepage  and  heat 
transfer  have  a  greater  effect  on  crack  length  than  does  friction 
alone  (for  the  current  parameters) .  This  effect  should  become  less 
dominant  at  later  times  since  there  is  an  inverse  time  relationship 
for  both  losses. 

5.4  Case  4  -  Specified  Borehole  Pressure  With  Friction  In  The  Fracture 
Included  (  e  =  4/xm  ) 

This  case  is  identical  to  Case  3  except  a  smaller  roughness  is 
used  in  the  friction  factor  calculations.  The  displaced  mesh  for  Case 
4  is  shown  in  Figure  5.41.  The  crack  is  about  7.3  meters  long  and  the 
crack  opening  diplacement  at  the  borehole  is  about  .0134  meters. 
Crack  length  and  opening  results  are  shown  in  Figure  5.42  and  5.43. 
KTn  during  propagation  is  shown  in  Figure  5.44.   Crack  tip  velocity 

results  are  shown  in  Figure  5.45.  Referring  to  Eqn  2.36,  it  can  be 
seen  that  a  factor  of  five  decrease  could  result  in  a  greater  than  two 
decrease  in  the  friction  factor,  for  predominantly  turbulent  flow. 
The  turbulent  term  should  dominate  for  the  cases  examined  in  this 
thesis.  The  decreased  roughness  has  a  minimal  effect  on  the  pressure, 
but  has  less  than  minimal  effect  on  the  temperature  when  compared  to 
Case  3.  Figures  5.46  through  5.50  show  the  gas  properties  at  various 
times  for  Case  4.  This  effect  is  due  to  the  fact  that  an  ideal  gas  is 
assumed  and  a  small  change  in  pressure  will  have  a  large  effect  on  the 
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Figure  5.36:  Gas  Pressure  for  Case  3  at  6.5,  9.5,  and  12  msec 
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Figure   5.37:    Gas  Velocity   for   Case   3   at   6.5,    9.5,    and  12   msec 
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Figure   5.38:   Gas  Temperature   for  Case  3  at  6.5,    9.5,   and  12  msec 
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Figure  5.39:  Gas  Density  for  Case  3  at  6.5,  9.5,  and  12  mse< 
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Figure   5.40:    Mach  Number   for  Case    3   at   6 . 5 ,    9.5,    and  12  msec 
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Figure  5.41:  Displaced  Mesh  Results  for  Case  4  at  12  msec 
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Figure    5.42:    Crack   Length   for   Case   4 
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Figure   5.43:    Crack  Opening  Displacement   for  Case   4   at   12  msec 
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Figure  5.44:  KID  During  Propagation  for  Case  4 
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Figure  5.45:  Crack  Velocity  for  Case  4 
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Figure   5.46:    Gas   Pressure   for  Case  4  at   6.5,    9.5,    and  12  msec 
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Figure  5.47:  Gas  Velocity  for  Case  4  at  6.5,  9.5,  and  12  msec 
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Figure  5.48:  Gas  Temperature  for  Case  4  at  6.5,  9.5,  and  12  msec 
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Figure  5.49:  Gas  Density  for  Case  4  at  6.5,  9.5,  and  12  msec 
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Figure  5.50:  Mach  Number  for  Case  4  at  6.5,  9.5,  and  12  msec 
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temperature.  Since  the  pressures  away  from  the  crack  mouth  are 
slightly  higher  for  Case  4,  the  temperatures  are  higher.  This  stored 
energy  will  have  an  effect  later  in  the  analysis.  However,  at  this 
point  in  time  the  similar  pressure  distributions  for  both  cases  will 
result  in  similar  crack  lengths.  The  crack  entrance  energy  plot  is 
shown  in  Figure  5.51. 

5.5  Case  5  -  Burn  Model  With  Friction.  Heat  Transfer.  And  Seepage 
Included 

Case  5  is  similar  to  Case  2  since  all  losses  are  included. 
However,  this  case  attemps  to  model  a  real  problem  by  using  a 
propellant  burn  model.  Propellant  properties  are  found  in  Table  5.2. 
The  displaced  mesh  is  shown  in  Figure  5.52.  At  12  msec,  the  crack  is 
about  4.0  meters  long  and  the  crack  opening  displacement  at  the 
borehole  is  about  0.0045  meters.  Plots  of  crack  opening  and  length 
are  shown  in  Figures  5.53  and  5.54.  The  borehole  pressure  history  is 
shown  in  Figure  5.55.  While  the  borehole  pressure  is  increasing,  the 
crack   tip  velocity,   Figure   5.56,   and  Kjn,   Figure  5.57,  are  also 

increasing.  After  the  borehole  pressure  peaks  at  60  msec,  the  crack 
tip  velocity  and  KTn  also  peak  and  decrease.   If  the  borehole  pressure 

were  allowed  to  decrease  further,  crack  arrest  would  occur.  To 
prevent  this,  the  borehole  pressure  should  probably  be  maintained 
above  50  MPa . 

The  gas  plots  are  shown  in  Figures  5.58  through  5.64.  As  is  typical 
of  the  cases  with  losses  included,  there  is  large  variation  of 
properties  from  the  crack  mouth  to  the  crack  tip. 
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Figure  5.51:  Crack  Entrance  Energy  for  Case  4 
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Figure  5.52:  Displaced  Mesh  Results  for  Case  5  at  12  msec 
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Figure   5.53:    Crack  Opening  Displacement   for  Case   5   at   12  msec 
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Figure  5.54:  Crack  Length  for  Case  5 
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Figure  5.55:  Borehole  Pressure  History 
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Figure  5.56:  Crack  Velocity  for  Case  5 
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Figure  5.57:  KID  During  Propagation  for  Case  5 
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Figure  5.58:  Gas  Pressure  for  Case  5  at  6.5,  9.5,  and  12  msec 
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Figure  5.59:  Heat  Flux  for  Case  5  at  6.5,  9.5,  and  12  msec 


Figure  5.60:  Mach  Number  for  Case  5  at  6.5,  9.5,  and  12  msec 
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Figure  5.61:  Gas  Density  for  Case  5  at  6.5,  9.5,  and  12  msec 
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Figure  5.62:  Gas  Seepage  for  Case  5  at  6.5,  9.5,  and  12  msec 
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Figure   5.63:    Gas  Velocity  for  Case   5   at   6.5,    9.5,    and   12   msec 
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Figure   5.64:    Gas  Temperature   for  Case   5   at   6.5,    9.5,    and  12   msec 
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Figure  5.65  shows  the  crack  entrance  energy  plot  for  Case  5  at  12 
msec.  Since  its  slope  is  less  than  that  for  Case  2,  it  seems  that 
less  of  the  available  energy  is  converted  to  kinetic  energy  than  is 
the  case  for  Case  2.  For  about  the  same  crack  length  as  found  in  Case 
2,  about  26%  less  energy  was  required. 

For  Case  5,  the  maximum  crack  tip  velocity  was  approximately  500 
m/sec .  To  have  a  successful  stimulation  of  100  meters,  the  crack  will 
need  to  be  driven  at  least  200  msec.  Figure  5.55,  in  addition  to 
experimental  data  (Schmidt  et  al,  1981),  illustrates  that  this  200 
msec  driving  time  will  not  be  the  result  of  using  conventional 
propellants  alone. 
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Figure  5.65:  Crack  Entrance  Energy  for  Case  5 
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CHAPTER  6 
CONCLUSIONS  AND  RECOMMENDATIONS 

6 . 1  Conclusions 

Cases  1  through  4  allowed  us  to  compare  the  effects  of  changing 
parameters.  The  most  striking  result  is  the  large  effect  of  friction, 
seepage,  and  heat  transfer.  Without  these  losses,  the  gas  flow  into 
the  crack  drives  the  crack  at  the  limiting  speed  of  the  rock  and  would 
lead  to  bifurcation.  When  the  losses  are  included,  gas  flow  into  the 
crack  is  reduced.  This  reduced  flow  limits  the  crack  velocity  to  the 
range  of  500-600  m/sec.  From  this,  we  can  estimate  the  necessary 
pressure  loading  time  to  drive  a  fracture  a  given  length. 

The  analysis  using  a  propellant  (Case  5)  allows  comparison  with 
field  tests  and  the  work  of  Nilson  and  Hanson.  The  most  significant 
result  is  that  the  predicted  crack  length  is  approximately  4  meters. 
This  is  much  smaller  than  the  100-200  meter  length  envisioned  for 
successful  stimulation.  This  fracture  length  is  consistent  with 
experimentally  observed  lengths.  As  was  mentioned  earlier,  no  field 
test  has  resulted  in  a  fracture  longer  than  10  meters.  Our  results 
from  Case  5  agree  with  this  experimental  result.  Nilson  (et  al., 
1985)  agrees  with  this  result  to  a  large  degree. 

These  calculations  and  supporting  experiments  show  that  long 
fractures  cannot  be  obtained  by  simple  propellant  burns.  The  energy 
stored  in  the  propellant  is  only  sufficient  to  drive  short  fractures. 
As   illustrated   in  Case  5,  the  crack  velocity  and  KT_  both  started  to 

decrease  after  the  propellant  was  past  its  peak  pressure.   To  maintain 
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crack  growth,  the  borehole  pressure  would  need  to  be  sustained  and  not 
allowed  to  drop. 

One  possible  solution  would  be  to  use  additional  propellant  (or 
other  stored  energy  source)  above  the  region  being  stimulated.  This 
would  allow  sufficient  energy  to  be  supplied  and  would  maintain  the 
borehole  pressure  at  a  high  level. 

6 . 2  Recommendations  For  Future  Work 

Since  the  pressure  loading  on  the  crack  is  the  dominate  feature 
of  the  tailored  pulse  stimulation  technique,  more  work  should  be 
applied  to  the  gas  model  solution  algorithm.  With  the  current 
method  (Lax) ,  it  was  necessary  to  apply  damping  and  a  time  reduction 
factor.  To  avoid  this,  other  solution  methods  such  as  upwind 
differencing,  Lax-Wendrof f ,  and  staggered  leapfrog  should  be 
considered.   Finite  element  implementation  is  also  a  possibility. 

To  model  a  complete  problem,  a  full  model  geometry  with  the 
borehole  in  the  center  should  be  used.  This  would  allow  the  use  of 
unsymmetrical  in- situ  stresses,  and  therefore  allow  crack  curvature  to 
occur.  With  this  type  of  model  geometry,  a  complete  multiple  fracture 
analysis  can  be  accomplished. 

In  addition  to  the  more  complete  model  geometry,  more  parameter 
studies  should  be  conducted.   The  effects  of  different: 

1.  in-situ  stresses, 

2.  permeabilities, 

3.  and  propellants 
should  be  considered. 
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ABSTRACT 

Gas-driven  fractures  are  of  practical  importance  in  the 
stimulation  of  gas  and  oil  wells,  rock  blasting  with  explosives,  and 
the  containment  of  nuclear  explosions.  There  is  at  present  no  general 
analytical  tool  to  predict  the  fracture  length  and  propagation 
direction  under  transient  loading  conditions. 

We  develop  a  coupled  model  of  gas -driven  dynamic  fracture  and  use 
the  model  to  analyze  well  stimulation.  The  model  consists  of  two 
parts:  (1)  a  finite  element  model  used  to  calculate  dynamic  crack 
propagation,  and  (2)  a  finite  difference  model  of  compressible, 
dynamic  gas  flow.  The  finite  element  model  assumes  that  the  fractures 
are  discrete  parts  of  the  mesh,  with  automatic  remeshing  performed  as 
the  cracks  propagate.  The  compressible  gas  model  is  developed  using 
the  conservation  of  mass,  energy,  and  momentum  equations.  Losses  due 
to  friction,  seepage  of  the  gas  into  the  crack  walls,  and  heat 
transfer  to  the  walls  are  included  in  the  analysis. 

Results  clearly  show  the  importance  of  the  loss  terms.  With 
losses  included,  the  crack  speed  is  limited  by  gas  flow  into  the 
fracture.  Calculations  using  a  propellant  burn  model  give  a  predicted 
crack  length  of  4  meters,  which  correlates  with  observed  field  tests. 


